RNA-Seq - differential expression using DESeq2

---

The Snf2 dataset

The RNA-Seq dataset we will use in this practical has been produced by Gierliński et al, 2015) and (Schurch et al, 2016)).

The dataset is composed of 48 samples of yeast wild-type (WT) strain, and 48 samples of Snf2 knock-out mutant cell line. The prepared RNA-Seq libraries (unstranded) were pooled and sequenced on seven lanes of a single flow-cell on an Illumina HiSeq 2000 resulting in a total of 1 billion 50-bp single-end reads across the 96 samples. RNA-Seq reads have been cleaned, mapped and counted to generated a count data matrix containing \(7126\) rows/genes.

In this tutorial, we will start from the count table generated by the authors of the study, and use different R tools in order to

• explore the distribution of the data;
• test alternative normalisation procedures;
• select differentially expressed genes (DEG);
• cluster the count table to regroup genes having similar transcription profiles;
• compare the lists of co-regulated genes with the annotations from the Gene Ontology, in order to test if they are statistically enriched for particular functional classes.

---

Loading the dataset

R enables to download data directly from the Web with the download.file function. After, we begin with some verification steps.

The expression matrix and phenotypic information will be loaded into R using the read.table function. Both table will be converted into a data.frame object when loaded into R.

• The countTable object should contain counts for each gene (row) and each sample (column).

• The phenoTable object should contain the designation of each genotype, that is 48 biological replicates of the wild type strain (WT) and the Snf2 mutant.

Note that the created tables include column names as factor object.

Reminder: a factor is a vector with levels (categories), which permits an efficient storage and indexing, but can in some cases lead to misleading effects. To circumvent this, we will sometimes have to convert the factor to a vector, with the R command as.vector().

# List of url and file names
url <- "http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data"
countFile <- "counts.txt"
phenoFile <- "expDesign.txt"

## Local paths: create a directory to store the data
dataDir <- ("~/jgb71e_practicals/rnaseq_snf2_Schurch_2015/data")
dir.create(dataDir, showWarnings = FALSE, recursive = TRUE)

# Download the data files unless they are already here
setwd(dataDir)

if(!file.exists(countFile))
  download.file(url = file.path(url,countFile), destfile = countFile )

if(!file.exists(phenoFile))
  download.file(url = file.path(url,phenoFile), destfile = phenoFile )

# Load the count table
countTable <- read.table(countFile, header=TRUE, row.names=1)
phenoTable <- read.table(phenoFile, header=TRUE, row.names=1)

# View part of the data/table
countTable[101:108, 1:5]

            WT1   WT2   WT3   WT4   WT5
yal034w-a    78   116   103   113   113
yal035w    2815  4160  2969  4411  4836
yal036c     331   599   348   525   691
yal037c-a     0    37    23    15    51
yal037c-b    70   183   199   208   196
yal037w      22    69    32    58    51
yal038w   21080 18869 14649 28551 20208
yal039c    1541  1447  1318  2251  1501

# Define a strain-specific color for each sample, and add it as a supplementary column
strainColor <- c("WT"="green","Snf2"="orange") # Choose one color per strain
## Define the column of each sample according to its genotype
phenoTable$color[phenoTable$strain == "Snf"] <- strainColor["Snf2"]
phenoTable$color[phenoTable$strain == "WT"] <- strainColor["WT"]

# Check the dimensions of the pheno table
colnames(phenoTable)

[1] "strain" "color" 

Exercice

1. Find the commands to visualize the first and last line of a data.frame object (clue: similar to bash commands)

View solution|Hide solution

head(phenoTable)
tail(phenoTable)

2. Check the completeness (dimensions) of the data (tables)

View solution|Hide solution

## Get the number of genes and sample from the dimensions of the count table
nb.samples <- ncol(countTable)
print(nb.samples)
nb_genes <- nrow(countTable)
print(nb_genes)

## Check that the dimension of the pheno table is consistent with the number of samples
dim(phenoTable)

3. Check the number of replicates in each class using the table function. To better understand what you could do with this function, play with it (varying the parameters) with the following test table.

## Quick test with a small table generated on the flight
testTable <- data.frame(
  ABC = sample(c("A", "B", "C", "D", "E"), replace=TRUE, 20),
  gender = sample(c("boy", "girl", "robot"), replace=TRUE, 20),
  yesno = sample(c("yes", "no"), replace=TRUE, 20))

View solution|Hide solution

## Check the content of the testTable
testTable

## Get stats about column content
table(testTable)
table(testTable[,2])
table(testTable[,2],testTable[,3])

## The same with column names instead of column numbers
table(testTable$gender)
table(testTable[, c("gender", "yesno")])

## Check the number of samples per strain in the yeast count table
table(phenoTable$strain)
table(phenoTable)

4. Draw a barplot showing the total number of reads in each sample (clue: use colSums function) with sample colored (according to the color column in phenoTable). What can you say from this diagram?

View solution|Hide solution

barplot(colSums(countTable)/1000000, 
        main="Total number of reads per sample (million)",
        col=phenoTable$color, 
        #        names.arg = "", 
        las=1,  horiz=TRUE,
        ylab="Samples", cex.names=0.5,
        xlab="Million counts")

Million reads per sample.

The barplot indicates the library sizes (total number of reads) for each sample. We can see important differences, ranging from 5.57 to 15.2 million reads per sample.

This has important consequences: any read count should be interpreted relative to the sequencing depth of the corresponding sample. We will thus need to normalise the counts. Consequently, it makes not much sense to compute summary statistics per gene (mean counts, standard deviation, …) before having normalized the data.

---

Descriptive statistics

Basic statistics

Before going further in the analysis, we will compute some descriptive statistics on the dataset. At this stage we only compute statistics per sample, since statistics per gene are meaningful only after library-wise normalization of the counts.

# Min, Max, median (...). 
# Here on the first 4 samples
summary(countTable[,1:4])

Among the reference functions in R, apply (and its derivates) allow to repeated execution of a function to a serie of data, at a noteworthy higher speed tha using for loops. Note that, if working on a matrix/data frame, the apply function requires that you indicate if you want to apply the function to the rows (MARGIN=1) or to the columns (MARGIN=2). Some built-in functions already exist as colMeans, colSums, rowMeans and rowSums. Let’s use apply to get the summary for every column of the count table. It produces a matrix as output that we transpose to obtain samples as rows and the statistics as columns, and move back to a data.frame object.

statsPerSample <- data.frame(t(apply(countTable, 2, summary)))
head(statsPerSample)

We can now add some columns to the stats per sample. For example, here are the library sums and the 5% quantile.

statsPerSample$libsum <- colSums(countTable) # library sums
statsPerSample$perc05 <- apply(countTable, 2, quantile, probs=0.05)

Exercice

Integrate some more columns to statsPerSample: - the 10%, 90% and 95% quantile;

View solution|Hide solution

statsPerSample[c("perc10", "perc90", "perc95")] <- 
  apply(countTable, 2, quantile, c(0.10, 0.90, 0.95))

• the number of zeros for each column;

View solution|Hide solution

statsPerSample$zeros <- apply(countTable==0, 2, sum)

• the percentage these zeros represent. View solution|Hide solution

statsPerSample$percZeros <- 100*statsPerSample$zeros/nrow(countTable)

Finally, randomly choose (with the sample() function) 10 samples (rows) of the stats per sample table and look at their statistics.

statsPerSample[sample(1:ncol(countTable), size = 10),]

## We can nicely format the result for the report, using the knitr::kable function
kable(statsPerSample[sample(1:ncol(countTable), size = 10),])

Distributions

The summary only displays a few milestone values (mean, median, quartiles). In order to get a better intuition of the data, we can draw an histogram of all values.

hist(as.matrix(countTable), col="blue", border="white")

Histogram of raw counts. The scale is determined by the gene with the highest count, which is apparently an outlier.

The histogram is not very informative so far, apparently due to the presence of a few very high count values, that impose a very large scale on the \(X\) axis. We can use a logarithmic transformation to improve the readability. Note that we will add pseudo counts to avoid problems with “zero” counts observed for some genes in some samples.

epsilon <- 1 # pseudo-count to avoid problems with log(0)
hist(as.matrix(log2(countTable + epsilon)), breaks=100, col="blue", border="white",
     main="Histogram of log2-transformed counts per gene", xlab="log2(counts)", ylab="Number of genes")

To get better insights into the distribution per sample, the boxplot offer a good perspective.

boxplot(log2(countTable + epsilon), col=phenoTable$color, pch=".", 
        horizontal=TRUE, cex.axis=0.5,
        las=1, ylab="Samples", xlab="log2(Counts +1)")

Box plots of non-normalized log2(counts) per sample.

Another way to get an intuition of the value distributions is to use the plotDensity function, which draws one density curve for each sample.

# We will require one function from the affy package
if(!require("affy")){
  source("http://bioconductor.org/biocLite.R")
  biocLite("affy")  
}
library(affy)
plotDensity(log2(countTable + epsilon), lty=1, col=phenoTable$color, lwd=2)
grid()
legend("topright", legend=names(strainColor), col=strainColor, lwd=2)

Densities of log2(counts). Each curve corresponds to one sample.

Beware: the R function plotDensity does not display the actual distribution of your values, but a polynomial fit. The representation thus generally looks smoother than the actual data. It is important to realize that, in some particular cases, the fit can lead to extrapolated values which can be misleading.

Scatter plots

Let’s have a look at the scatter plots using the pairs function. We will only represent 10 randomly selected samples.

# Define a function to draw a scatter plot for a pair of variables (samples) with density colors
plotFun <- function(x,y){ 
  dns <- densCols(x,y); 
  points(x,y, col=dns, pch=".", panel.first=grid());  
#  abline(a=0, b=1, col="brown")
  }

# Plot the scatter plot for a few pairs of variables selected at random
set.seed(123) # forces the random number generator to produce fixed results
pairs(log2(countTable[,sample(ncol(countTable),5)] + epsilon), 
      panel=plotFun, lower.panel = NULL)

Scatter plot of log2-counts for a random selection of samples.

The command pairs draws a scatter plot for each pair of columns of the input dataset. The plot shows a fairly good reproducibility between samples of the same type (WT or KO, respectively): all points are aligned along te diagonal, with a relatively wider dispersion at the bottom, corresponding to small number fluctuations.

In contrast, on all the plots comparing a WT and a KO sample, we can see some points (genes) discarding from the diagonal.

---